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ABSTRACT 

A multigrid method is presented for calculating turbulent jets in crossflow. Fairly 
rapid convergence is obtained with the k — e turbulence model, but computations 
with a full Reynolds stress turbulence model (RSM) are not yet very efficient. Grid 
dependency tests show that there are slight differences between results obtained on 
the two finest grid levels. Computations using the RSM are significantly different from 
those with k - e model and compare better to experimental data. Some work is still 
required to improve the efficiency of the computations with the RSM. 
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constant in Reynolds stress model " “ ~ 
k — e turbulence model constants 
jet diameter 

near-wall proximity function in Reynolds stress model 
rate of production of turbulent kinetic energy 
height of duct 
turbulent kinetic energy 
pressure 

jet to crossflow velocity ratio 
pitch spacing of multiple jets 
source term for dependent variable $ 
crossflow velocity 
cartesian velocity components 

Reynolds normal stresses in cartesian directions 

Reynolds shear stresses 

jet velocity 

width of duct 

cartesian coordinates 

constants in Reynolds stress model 

Kronecker delta 

rate of dissipation of turbulent kinetic energy 

von Karman constant 

thermal/species diffusivity 

molecular viscosity 

turbulent eddy viscosity 

density 

turbulent Prandtl/Schmidt number for $ 

General representation of dependent variable 


Superscripts 

1 lateral direction 

2 vertical direction 

3 longitudinal direction 
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INTRODUCTION 


Three-dimensional turbulent jets in crossflow have important engineering applications in 
both confined and unconfined environments. Examples of jets issuing into confined cross- 
flow include internal cooling of turbine blades, dilution air jets in combustion chambers, 
jets from V/STOL aircraft in transition flight, etc. The examples of turbulent jets issuing 
into unconfined (semi-infinite) crossflow are even more numerous. These include discharges 
from cooling towers or chimney stacks into the atmosphere or sewerage and waste heat 
into water bodies, film-cooling of turbine blades, etc. 

The interaction of the jets with the crossflow has been investigated in numerous exper- 
imental studies [1-4]. Crabb et al [2] present a comprehensive review of pre 1980 studies, 
most of which only deal with mean flow properties. Measurements of turbulent properties 
can be found in [2-6]. Numerous computational studies of the generic problem of turbulent 
jets in crossflow are also reported in the literature [7—10], but there are uncertainties as to 
the accuracy of the results. First, it has not been possible to obtain grid independent solu- 
tions. In a recent study, Claus and Vanka [11] present results with computational grids up 
to 96x96x256, but could still not confirm grid independency. Secondly, most computations 
have employed the Jfc — e turbulence model which assumes gradient diffusion relations for the 
Reynolds stresses and an isotropic eddy-viscosity distribution. Measurements of Reynolds 
stresses and velocity fields by Andreopoulos and Rodi [4] suggest that these assumptions 
may not be appropriate in many cases. 

The present study tries to address these problems. Computations are performed with a 
multi-grid procedure which enables convergence on very fine grids within a relatively small 
number of iterations. A second-moment closure turbulence model is utilised to compute 
the Reynolds stresses, as well as the more popular k — e model. 


MATHEMATICAL MODEL 


In the present work we solve the time-averaged, three-dimensional, steady state equations 
governing the turbulent flow and heat transfer. The equations may be expressed in carte- 
sian tensor notation as: 
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with i=l,2,3 and j=l,2,3 representing properties in the lateral, vertical and longitudi- 
nal directions, respectively, y ’ (= y 1 , y 2 , y 3 ) represents the cartesian coordinates; 17, the 
cartesian velocity components; P the presssure and $ the normalised temperature or con- 
centration. p is the density, p is the molecular viscosity and A is the thermal or species 
diffusivity. The equations are expanded by using Einstein’s summation rule for repeated 
indices. -puTtZJand -puT^are respectively, the Reynolds stresses and heat/concentration 
fluxes which must be determined by a turbulence model before the system of equations 
can be closed. 

Turbulence Model 

In this paper results are presented of calculations with both the standard k — e turbulence 
model [12] and a second-moment turbulence closure based on the proposals of Launder, 
Reece and Rodi [13], hereafter denoted LRR. 

In the standard k — e model, the Reynolds stresses and heat fluxes are approximated 
with the Boussinesq eddy viscosity /diffusivity concept as 



, dUj\ 2.C 
+ Ihf ) 3 k6ij 

( 4 ) 

pui<p = 

dy % 

( 5 ) 


/i t is the eddy viscosity given by : 

Ht = ( 6 ) 


The distributions of the turbulent kinetic energy k and its rate of dissipation e are then 
obtained from the solution of the transport equations: 

W {pU ’ k) = W (ft#) +G - ( ’ 6 < 7 > 

^7 ( pU i e ) = ^ + C “ G ! “ M 

where G is the turbulence production rate given by: 

G = -puTtZJ (9) 

The empirical constants appearing in the equations above are: 
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( 10 ) 


c M = 0.09, c el = 1.44, c £2 = 1.92, a t = 0.9, <7* = 1.0, cr £ = 1.3 


and Sij is the Kronecker delta. Equations (l)— (9) form a closed set which is solved itera- 
tively in a sequential manner to obtain the velocity and temperature fields. 

For the second-moment closure model, we adopt the proposals of LRR [13], model 1 
to approximate the pressure-strain, diffusion, and dissipation terms in the Reynolds stress 
equations. The resulting system of equations can be written in cartesian tensor notation as: 
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with the empirical coefficients a, /?, 7, Ci, and c, given by: 


a = 0.7636-0.06/;/? = 0.1091 + 0.06 /; 7 = 0.182; c k = 1.5-0.50/; c 3 = 0.22 

( 12 ) 


and f is the wall-proximity function with value of unity near walls and zero in a completely 
free stream. 

To be consistent with the second-moment closure, we modify Eq. (8) slightly by intro- 
ducing non-isotropic eddy-viscosity distributions into the diffusion terms, following LRR. 
We do not apply a second-moment closure model to the turbulent heat flux u^p in the 
present study. Equations (1), (2), (3), (8) and (11) form a closed set which should be 
solved simultaneously. The equations are solved in a sequential manner, first for the ve- 
locity components and pressure, then Reynolds stresses and dissipation, and finally the 
temperature. If the terms involving gradients of the Reynolds stresses on the r.h.s. of Eq. 
(2) are treated explicitly the system of equations will be very stiff and it will be extremely 
difficult to obtain a converged solution with an iterative scheme. The stiffness can be 
reduced considerably by splitting the Reynolds stress u JuJ into two parts: 
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The first part is treated explicitly. The second part is added to the molecular diffusion 
term and treated implicitly. The modified momentum equation has the form: 
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The set of Eqs. (1), (3), (8) (11) and (14) is solved, using underrelaxation factors of 0.75 
for the three velocity components, 0.25 for the presssure, and 0.7 for the temperature and 
all turbulence quantities. Four types of boundary conditions are encountered, namely: 
inlet, outlet, symmetry and walls. Inlet conditions are specified from experimental data. 
The outlet is an outflow boundary requiring no formal specification of conditions. Along 
symmetry planes the normal gradients of all variables are set to zero, and the normal 
velocity component is also zero. The walls are special in that we do not integrate all the 
way down, rather we use the wall-function method [12,13] to prescribe the values of the 
dependent variables at near-wall nodes. 

Multigrid Procedure 

In the present work the FAS-FMG (full approximation storage-full multigrid) algorithm 
originally developed by Brandt [14] is employed to solve the hydrodynamic equations. 
The present implementation derives from previous works by Demuren [15] and Vanka [16]. 
There are however significant differences. First, the present method uses a regular grid sys- 
tem with no staggering of the velocity nodes relative to the pressure nodes. The expected 
odd-even decoupling problem is overcome by adding a fourth-order artificial dissipation 
term to the pressure gradient. It can be shown that, with a coefficient of unity, this 
practice is equivalent to the so-called “momentum interpolation” method of Peric [17]. 
However, there is now the flexibility to vary the coefficient all the way down to zero, if nec- 
essary. The second difference is that the system of equations is now solved in a sequential 
manner as opposed to the coupled approach proposed by Vanka. Numerical experiments 
showed no advantage in using the latter in a multigrid procedure, and it can be shown 
mathematically that it is less stable in a single-grid procedure. Further, it appears easier 
to vectorize the decoupled procedure. 

The smoothing scheme is based on the SIMPLEC method described by van Doormaal 
and Raithby [18]. The salient steps are: 

1. Solve the Ui momentum equation using a guessed pressure field. 

2. Then the Ui momentum equation. 

3. Then the U s momentum equation. 

4. Compute the velocities on the faces of the control volume, each by linear interpolation 
plus a fourth-order artificial dissipation term. 

5. Then compute the mass source error in each control volume. 

6. Solve a pressure-correction equation to eliminate the mass source errors, and then 
correct the pressures and corresponding velocity components. 
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7. If on the finest grid level, solve the equations for k, e, and u.u,, as the case may be. 

8. Then solve the temperature equation. 

These steps are repeated until convergence. A V-cycle multigrid algorithm is utilised 
with 10 iterations on the coarsest grid and 3 on intermediate grids. One or more iterations 
are carried out on the finest grid depending on the achieved smoothing rate. Restriction 
from fine to coaxse grid is by averaging, and prolongation from coarse to fine grid is by 
trilinear interpolation. An ADI scheme is employed to solve the final set of algebraic 
equations for all variables at all grid levels. The underlying algorithm is the tri-diagonal 
matrix algorithm (TDMA) which is known to be recursive, and would thus not normally 
be vectorizable. However, by changing the data structure we can make all the internal 
loops of the ADI solver vectorizable on the Cray computers. Although we cannot remove 
the recursivity of the algorithm, the change in data structure ensures that all floating point 
operations are in vector form. Typical saving in total CPU time resulting from this is of 
the order of 50%. 

The equations for turbulent quantities k, e and uiuj are solved only on the current finest 
grid during the iteration process. Corresponding operators on coarser grids are calculated 
using restricted values for these quantities. However, the solution process on any fine grid 
is started with variable values prolongated from the converged solution on the immediate 
coarser grid. This point is especially important for the Reynolds stresses. 

RESULTS AND DISCUSSION 
Computational Details 

The test cases for the present work are derived from experimental studies of opposed 
jets discharging normally into a cross-stream by Atkinson et al [3]. Figure 1 shows the 
schematic diagram. Three cases are considered; two with a single pair of opposed jets and 
one with a row of five pairs. The computational details are presented in table 1, below. 

Table 1 :Computational Details of Test Cases. 


Case 

R 

S/D or W/D 

Coarsest 

grid 

Total Points 
on Finest Grid 

CPU Time/ Work Units 
(Cray YMP mins) 

1 

1.8 

12 

12x10x22 

876,744 

39/82 

2 

1.0 

12 

12x10x22 

876,744 

34/71 

3 

1.8 

4 

8x10x22 

534,600 

27/94 


In table 1, R represents the jet to crossflow velocity ratio, D is the jet diameter, W 
is the channel width in the single-pair jet study, and S is the jet-centreline pitch spacing 
in the multiple-pairs jet study. The channel heigth is 4D in all cases, so that there is a 
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symmetry plane at 2D. Cases 1 and 2 also have a symmetry plane passing through the 
centre of the jet and channel, whereas case 3 has a third symmetry plane through the 
mid-plane between jets. The inlet plane for the crossflow is located 4D upstream of the jet 
inlet and the outlet plane 14D downstream. We take advantage of the various symmetry 
planes to minimize the size of the computational domain. Figure 2 shows typical velocity 
vectors in the center plane computed with the RSM model on a 3-level grid. Details of 
the flow field are discernible in the exploded view (a). We can see the complex flow spiral 
in the wake, which shows how the crossstream passing between the jets is entrained into 
the jet. There is also a small reverse flow near the stagnation point in front of the jet. 
Further, for this velocity ratio, it is clear that the two opposing jets impinge upon each 
other along the symmetry plane, from about one diameter downstream of the exit pipes. 
The performance values in table 1 are for convergence to a normalised residual level of 
5xl0 -4 . They represent the total computational work done on all grids for computations 
with the k — e turbulence model on a single processor. Computations with the Reynolds 
stress model typically require twice as many iterations. It now appears that the multigrid 
scheme should also be applied to the turbulence quantities k, e and u,uy. The computations 
presented in the paper are for 3- and 4-grid levels using the k — e turbulence model and 

3- grid levels using RSM. 3-grid-level calculations take roughly 3-4 minutes. 

Grid Dependency 

One of the main aims of the present study is to attempt to obtain grid-independent results 
through considerable grid refinement, so as to separate numerical errors from model errors. 
Figures 3-5 present the results of grid dependency tests. In Fig. 3, computed longitudinal 
velocity contours are compared along four axial planes for case 1. The computations on 
the right are on a grid which is exactly twice as fine as the one for the results on the left. 
We see that the flow patterns are exactly the same, but there are slight differences in peak 
values, especially in the near field. Figure 4, showing computed concentration/normalised 
temperature contours emphasizes these points. We also see steeper gradients near the edge 
of the jet for the finer grid computations, which is a result of reduced numerical diffusion. 
Figure 5 repeats these comparisons along two axial planes for case 3, which has the same 
velocity ratio as case 1 but for multiple pairs of jets. It is interesting to note that the 
results for these cases are similar, except for the higher velocity magnitudes and smaller 
jet spreading in case 3. These are of course due to the higher blockage factor in this case. 
Although the results cannot be claimed to be grid-independent the differences between 
them are small enough for us to estimate changes that could be expected from further grid 
refinement as compared with turbulence model changes. These results are for 3-level and 

4- level grids. The next level of grid refinement will contain about 6.8 million points and 
should require about 256 megawords memory and 4 hours computational time. Although 
this is still within the capacity of the computer, the turnaround is expected to be very 
long. 
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Comparison with Experiment 

The present computations on a 3-level grid are compared with experimental results of 
Atkinson et al [3] in Figs. 6 and 7. The choice of grid is dictated by the little differ- 
ences in predictions between 3- and 4-level grids shown in the preceeding subsection, and 
the much greater expense of computations with RSM. Figure 6 presents the longitudinal 
velocity results for case 1 and fig. 7 those for case 2. For case 1, predictions with the 
k — e model at y*/D = 4 show the characteristic kidney shape. This is not evident in 
the measurements, and the RSM results show a much milder shape. It appears that the 
k — c model produces counter-rotating vortices that are too strong, leading to excessive 
distortion of the velocity contours. The Reynolds stresses act to reduce this motion and 
thus procure better agreement with the measurements. This trend is repeated in the re- 
sults at y 3 /D = 8 .In general, the RSM model produces contour shapes that are much 
closer to the measurements than the k — e model. The magnitudes may not be exactly the 
same, but the trends are very encouraging, and the changes are certainly more significant 
than obtained through grid refinement. However, more work needs to be done to improve 
the RSM predictions. Concentration comparisons are not presented because these are not 
truly sy mm etrical in the experiments since only one jet was seeded in each case. 

CONCLUDING REMARKS 

A multigrid procedure for calculating turbulent jets in crossflow has been presented. The 
procedure is applied with either a k — e turbulence model or a second-moment closure 
model. With the k — e model version, convergence can be obtained in less than 100 
equivalent fine grid iterations on any grid. RSM version still needs some work to bring it 
to the same level of efficiency. Computations with the k — e model show slight changes 
with grid refinement, but the use of the RSM leads to much closer agreement of velocity 
contours with experimental data. 
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FIGURE 1. - MULTIPLE PAIRS OF OPPOSED JETS IN CR0SSF10W. 
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(b) WHOLE FIELD. 

FIGURE 2. - COMPUTED VELOCITY VECTORS IN CENTER-PLANE, R = 1.8, 
W/D = 12. 
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